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Abstract 

We consider a reaction-diffusion process with retardation. The particles, immersed in traps 
initially, remain inactive until another particle is annihilated spontaneously with a rate A at a 
certain point x. In that case the traps within a sphere of radius R(t) = vt a around x will be 
activated and a particle is released with a rate \x. Due to the competition between both reactions 
the system evolves three different time regimes. While in the initial time interval the diffusive 
process dominates the behavior of the system, there appears a transient regime, where the system 
shows a driveling wave solution which tends to a non-trivial stationary solution for v — > 0. In 
that regime one observes a very slow decay of the concentration. In the final long time regime 
a crossover to an exponentially decaying process is observed. In case of A = /i the concentration 
is a conserved quantity whereas for fi > A the total particle number tends to zero after a finite 
time. The mean square displacement offers an anomalous diffusive behavior where the dynamic 
exponent is determined by the exponent a. In one dimension the model can be solved exactly. In 
higher dimension we find approximative analytical results in very good agreement with numerical 
solutions. The situation could be applied for the development of a bacterial colony or a gene-pool. 
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I. INTRODUCTION 



A broad variety of different problems in physics and biology can be formulated in terms 
of generalized evolution equations including delay and feedback effects . Such memory- 
controlled effects should be a further unifying feature of complex physical 0, 0] as well as 
biological systems P far from equilibrium. Recently [f| memory effects have been qual- 
ified as the key feature in describing dynamical systems. Distributed delays are able to 
stabilize ecological systems 3). Even the traffic flow can be characterized by the scaling 
properties of the underlying memory function ^, |^. The whole history of systems offering 
self-organized criticality as earthquakes is discussed in Whereas most of the papers 

are addressed to a purely time dependent but homogeneous memory, the present one gives 
an extension to spatiotemporal processes. Based on our recent studies of several evolution 
models [ill . 12,Q,Q], where such non-Markovian memory effects had been enclosed, and 
on the analysis of chemical reactions with a long range memory H, the present study is 
focused to a diffusion-reaction behavior with short and long range memory couplings. As an 
essential new element we include the case that retardation effects are relevant. Thus cells 
are able to interact as well as by short-range forces such as adhesion and through either 
long-range forces or retarded interaction via chemical signals la ]. These effects have an 
essential influence on the long time behavior, but they have also an impact on the dynamics 
in an intermediate time regime. As the consequence of the interplay between conventional 
diffusive transport and annihilation and creation processes with feedback we demonstrate 
the system evolves slowly varying plateau in the concentration. 

Our analysis can be grouped into the continuing interest in feedback couplings. A memory 
dominated behavior is well established in analyzing the freezing processes in under-cooled 
liquids 3, 3> where the underlying mathematical representation is based on a projector 



formalism proposed by Mori [2fJ. As the result of the projection procedure the irrelevant 
degrees of freedom, which are rapid fluctuating ones, contribute to the time evolution of 
the relevant degrees of freedom as well as by instantaneous and by delay-controlled terms. 
Since the projection formalism is not restricted to selected systems, the modification of the 
evolution equation due to memory effects is rather generic. The relevant degrees of freedom 
offering a slowly varying behavior could be the concentration of certain species or the proba- 
bility density for the particles. The crucial influence of memory effects can be illustrated by 
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considering a single particle moving in a disordered environment. Due to the strong disor- 
der each member of an equally prepared ensemble makes experience by its own surrounding, 
which is even modified by the random walk of the single particle itself. Such a behavior can 
be described by a non-linear Fokker-Planck equation with an additional memory term as 
had been demonstrated in [21]. The analytical results, based on a one- loop renormalization 
group analysis, is in accordance with numerical simulations (h| and further analytical stud- 



ies 



12|. Whereas in that approach the memory effects are believed to be originated by the 



inherent non-linear interaction of the many particle system themselves, i.e. the time scale 



of the memory is determined by the relevant varia 



Die itself, there is a broad class of models 



22, 



23 



241 ]; for a survey and applications in 



which are subjected to external delay effects 
biology, see 21 

Recently |26( memory effects in correlated anisotropic diffusion are studied in nanoporous 
crystalline solids. Likewise the effects of transport memory are discussed within the Fisher's 
equation fl, ratehet devises Q also applicable for bacterial population dynamics Q. 
There appear a non-linear damping and traveling wave solut.ons fa* The transport with 
memory, depending on the survivability of a population, is analyzed in [3JJ . 
In case the transport is realized after a spatiotemporal accumulation process the time evo- 
lution of the probability or the concentration of a different species could be also dependent 
on the history of the sample to which it belongs. Thus the evolution of the relevant variable 
p(x, t) has to be supplemented by memory terms. One possible but rather general structure 
of an evolution equation with memory should read j^J 

t oo 

d t p(x,t) = M(x,t;p, Vp) + J ' dt' J d d x'JC{x-x',t-t';p,Vp)C(x',t';p,Vp). (1) 

-oo 

This equation is of convolution type. Whereas the operator M. characterizes the instanta- 
neous and local processes, both operators K and C are responsible for the delayed processes. 
In general the operators M., K, and £ may be also non-linear in p(x, t) and Vp(x, t). They 
have to be specified according to the physical situation in mind. In particular we have 
shown [2, 0] that the form of the operator C is restricted when p(x, t) is assumed to be 
conserved. Regarding the large variety of systems with feedback couplings it seems to be 
worth to study simple models, which still include the dynamical features of evolution models 
as conservation of the relevant quantity p(x, t) and moreover, a time-delayed coupling. 
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In the present paper we extend the analysis by considering another kind of memory not 
included in Eq. (JIJ . To that aim let us consider a diffusion- reaction model with retardation 
for a concentration p(x,t). Different to Eq. there appears a spatiotemporal coupling 
well known from retarded potentials. Each particle makes a random walk described by a 
diffusion equation (^tJ in the continuum limit. Moreover, it can be annihilated with the rate 
A spontaneously. As a competing process a self-organized creation of particles with rate /i 
is introduced. However the additional process is subjected to a retardation. The situation 
corresponds to a diffusion-reaction process where the particles are able to emit permanently 
signals. After a finite running time inactive particles will be excited. 



II. DIFFERENT REALIZATIONS AND THE GENERAL SOLUTION 

To be specific let us consider a diffusion-reaction model for the concentration p(x, t) of 
the form 

= DV 2 p(x,t) - Xp(x,t) 

+ fi 1™ d d x'<$>(x',t- \ x ~ x '\ \ N(\ x-x' |)6 (t- \ x ~ x ' H (2) 



Particles, performing a random walk with diffusion constant D are annihilated with the rate 
A at a certain point x at time t. The subsequent term describes the creation of particles with 
the rate \i. This creation process is triggered by a signal that will be transmitted by each 
particle at the time t' = t — JEzhLI. With other words, particles are created whenever the 
signal wave, originated from any point x' at time t' is reached the point with the coordinate 
x at the observation time t < t' . The quantity $>(x',t') is determined by the concentration 
p, and different realizations will be discussed below. To complete the model we have to fix 
the normalization factor N(r) by J r d ~ 1 N(r)dQ ( i = 1 which leads to 

I 2vr d / 2 
N ( R )=rT5tt With ^ 



n d R d -' a r(d/2) 

After Fourier transformation Eq. (j2J reads 
dp(k, t) 



dt 



+ [Dk 2 + X]p(k,t) = fiT(d/2)2 d / 2 - 1 [ dt'[kv{t-t'} l ' d/2 Jd/2-i{kv{t-t')^{k,t') 

Jo 
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where J a (x) is the Bessel function. Performing Laplace transformation with respect to time 
we get 

(z + Dk 2 + X)p(k, z) = p(k, t = 0) + ^F Q,l, t -^pj z) . (3) 

Here, F(a,b,c; y) is the hypergeometric function. For the further discussion the function 
t) has to be specified. Here we discuss three cases 

(i) $(£,t) = p(x,0)5(t)= Po (x)5(t) 

(a) $(£, t) = p(x,t) 

i(x 

In case (ii) the concentration is given by 



in) $(f,t) = ^M. (4) 



FV ' ; z + Dk 2 + X + fF(l,l,l;-!^f) ' 

In case (iii) the retardation is manifested by a coupling to the time variation of the concen- 
tration itself. It results 



Po(k) 



^ + DA; 2 + A + fF(|,l,f;-^) ' (6) 

Both equations do not give rise to an transient regime as it will be observed in the first case. 
Therefore we restrict the further discussion mainly to the first case which is discussed in 
detail in the forthcoming section. 



III. LONG RANGE MEMORY 



In case of 

$(x,t) = Po (x)5(t) (7) 

the model offers additional specifications like conservation and anomalous behavior. The 
realization given by the last equation means that the creation process is triggered by a signal 
that had been transmitted at the initial time t — 0, where the particles are distributed with 
the concentration po(x). The signal is presumed to be permanent existent in the whole space 
like the microwave background radiation in cosmology. The activation of particles by this 
background signal can be further generalized in the following manner. If at an arbitrary 
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point x a particle is annihilated then the activation process sets in. However, to simulate 
the influence of the backward signal it should be restricted to a sphere with radius 



R(t) = vt a (8) 

around the mentioned point at x, where the exponent will be specified below. Only in case 
a = 1 the quantity v can be interpreted as a real velocity. While the particles perform a 
random walk, they will be annihilated with a certain rate A and created with the rate p. 
This process is governed by signals coming from the a sphere with radius R(t) around the 
point with the coordinate x. The situation in mind is simply described by the following 
evolution equation 



Of 



DX7 2 p(x, t) - J d d r [ \p(x, t) - p N(r)p( x - f, 0)] 5{r - R[t)) with r =\r\ . 

(9) 

Here the first term describes conventional diffusion with the diffusion constant D. The 
second term stands for the annihilation of particles with rate A. The last terms models the 
situation discussed before. Particles are created with a rate p whenever the signal wave 
originated from the initial state of the system is reached the spatial point x. This process 
corresponds to a kind of self-stabilizing of the system, especially for A = p, where the total 
particle number is conserved. If both rates are different one obtains for the total number of 
particles P(t) the relation 

Pit) = P(0) + e- Xt ( 1 - 1 , (10) 



-A V A, 

which is also independent of the realization of the function R(t). In the limiting case A = fi 
the conservation is guaranteed. In case of a vanishing rate A = there occurs a linear 
increase of the concentration according to P(t) = P(0)(1 + fit). Let us remark that in the 
opposite case where growth and death rate will be interchanged, i.e. A — > —A and p, — ► — p 
the last equation leads to a complete extinction of the species after a finite time 

< X ln^ " 



A \p — A 

provided p > A. As one can see from Eq. © there appears a driving force between the 
concentration of the spatial point x and a point far from that at x—R(t), i.e. all spatial points 
immersed in a sphere of radius R(t) = vt a contribute to the instantaneous development of the 
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system. The model describes a situation, where a set of particles performing a random walk, 
carry the concentration in the past like a jet. In case of zero parameter v — the system 
reveals an extreme long range memory, characterized by a direct coupling of the probability 
p(x,t) to its initial value p(x, 0). Due to this long range coupling there exits a non-trivial 
stationary solution as demonstrated in This stationary solution is originated by the 

balance of the diffusive current and the memory induced current proportional p(x, t)—p(x, 0). 
In the present approach we consider a new situation consisting of the existence of an active 
local environment of radius R(t) around each spatial point at x. The situation corresponds 
to the presence of traps in which particles are immersed. These traps will be installed 
at initial time t = 0. However the barrier to overcome is to high, so that the particles 
are confined within the traps. The particle inside a trap remains inactive until another 
particle is annihilated spontaneously at a spatial point x. Let us assume that through the 
annihilation process energy will be released so that the particles are able to overcome the 
barrier of the traps and the particle included become active with a rate /i provided they are 
in a sphere of radius R(t) around x. With other word, instead of a long range memory there 
is a short range one induced by a coupling of the current concentration p(x, t) to the value 
of the distribution p(x — R(t), 0). In that case we demonstrate that the stationary solution 
disappears and instead of that the system evolves a very slow decay in a broad transient 
time regime. In the limiting case of vanishing v the stationary solution is restored. The 
transient regime is likewise originated due to the balance of the diffusive current and the 
additional memory current proportional to p{x,t) — p(x — R{t)). Assuming the quantity 
p(x, t) is proportional to the probability to find the particle at the spatial point x at time t 
we get for the mean square displacement s(t) = (x 2 ) the evolution equation 

= 2dD - Xs(t) + n[s(0) + R\t)} . (11) 

In case of R(t) = vt a the last equation can be solved resulting in 

s( t ) = ^ (1 - e- xt ) + fiv 2 (-X) 1+2a 7 (1 + 2a, -Xt) (12) 
A 

Here 7(0, y) is the incomplete gamma function. For nonzero A the asymptotic behavior of 
the displacement is for At ^> 1 given by 

, , t 2a 
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while in the case of vanishing annihilation rate A = one gets exactly 

2 

lim s(t) = 2dDt + Jf^t 1+2a . 
a^o w 1 + 2a 

The growth term in Eq. (J2J) gives rise to anomalous diffusion which dominates for a > the 
long time behavior of the system. In particular, even for a = 1/2 globally the system offers 
a ballistic behavior. In case of conservation of p(x, t) the last relation leads to diffusion due 
to A = \i. 

IV. SOLUTION 

Since the system exhibits no general analytical solution for an arbitrary behavior of R(t) 
we restrict the further discussion to the case a — 1. Fourier transformation of Eq. (j2J yields 

_ Jd=i(kR(t)) 

d tP (k,t) = -(Dk 2 + \)p(k,t) + XA Po (k)^ ^ with A = T(d/2)2 d l 2 - 1 - k=\k\ . 

(13) 

Here J„(y ) is the Bessel function. Although the last equation is only valid for dimension 
d > 1 |3Ja], a separate transformation reveals that the one dimensional case with 

d tP (k, t) = -(Dk 2 + X)p(k, t) + fipo{k) cos(kR(t)) , (14) 

is also included. To analyze whether the model offers a nontrivial stationary solution it 
is appropriate to consider the Laplace transformation p(k,z). Aside of special cases the 
transformation can be fulfilled only for R(t) = vt, i.e. for a — 1. From Eq. (jlHj) we find 

where F(a, b, c; y) is the hypergeometric function. The result is always valid for d > 1. The 
long time limit results for z — > 0. Because there is no general asymptotic behavior of the 
hypergeometric function F(a, b, c, y) for y —>■ oo and arbitrary dimension d we will discuss 
different dimensions separately. 
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A. Exact solution for d = 1 



In the one- dimensional case the solution for p(x, t) can be expressed by elementary 
functions. In d = 1 the evolution equation reads 

dp(x,t) d 2 p(x,t) 



D- 



()f - -\p(x,t) + £\p(x-R(t),0)-p(x + R(t),0)] with p(x,0)= Po (x), 

(16) 

where R(t) = vt is defined by Eq. (JHJ) with the exponent a — 1. Performing Fourier 
transformation with respect to the spatial coordinate we can solve the the resulting equation 
for p(k,t). According to Eq. (JT3J) we get 



p(k,t) = po(k)e 



{Dh?+\)t + Po( k )^ ( £ 



ikvt _ -(Dk 2 +X)t 



Dk 2 + ikv + X 



+ c.c. 



(17) 



By solving this equation one observes that the three terms describe three different time 
regimes. Whereas the first one, as solution of the homogeneous equation, represents the 
initial time interval, the second term gives rise to a traveling wave behavior which tends 
for v = to a space dependent stationary solution. The last term is responsible for the 
final time interval and leads to an exponential decay. After a lengthy but straightforward 
backward transformation we find p(x, t) = Pi(x, t) + p s (x, t) +pf(x, t) with 
„-xt 



Pi(x,t) 



p s {x,t) 



p f (x,t) 



dx po(x')e 4 Dt 



2(JL 



dx'p (x') \e(x -x' + vt)e~ Ul{x - x ' +vt) + Q(x' -x- vt)e~ U2{x J ~ x ~ v€) 



fie' 



-xt 



n 



dx'po(x') I e-M—'-Dtu, 



)erfc 



+ e 



"2 



(x-x>+Dtv 2 ) eTfc f 2v 2 Dt + (x - x') \ 1 

V vwt J S 



with Q = vV + AXD, 



Vl/2 



n =f v 
2D 



> 0. 



(18) 



In the limiting case v = the part p s (x,t) tends to a non-trivial stationary solution 



Ps{X) 



'\D 



+ 00 



| x — x 1 1 



' D 



dx'p (x')e « with ^ = \j — . 

A 
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This result is in accordance with a previous study 17[. The behavior of p(x,t) is controlled 
by a single dimensionless parameter 

v 



VXD 7 

which is the ratio of the two velocities v and V XD. 
B. Solution in higher dimensions 

Guided by the solution in d = 1 let us analyze the three dimensional case. In d = 3 the 
hypergeometric function in Eq. (fTKj) can be expressed by an elementary function resulting 
in 

,r* s ,? A + v/kv arctanf— ) , 

Although we are not able to get a general solution for p(x, t), one can estimate the solution 
in different regimes. Denote I a typical length scale and r a characteristic time scale we can 
distinguish two cases: (i) v ^> - and (ii) v -C -■ Due to — 3> 1 in the first case the last 
equation leads to 

p(k,z)*p (k) *tflx - (2°) 
z + Dk l + A 

From here we get an exponential decay p(k, t) oc exp(— [Dk 2 +X]t) . In more detail we obtain 

-At 

' j3 / 



£) « ^ + p f = ^ nDt j 3/2 J d xp (x)e 

+ £^ / ^WK^F, (i/2,3/2; <*z|>!) . (21) 

Here i-Fi is the degenerate hypergeometric function which can be expressed by elementary 
functions: X F X (1/2, 3/2, y) = ^/tt/Av erfi( v ^) Q- Notice that due to the condition If > 1 
the result can be not extrapolated to v = 0. In the opposite case — <C 1 Eq. (|19p is rewritten 
in the following form 

p(jfe,z) « ^® + tf(jfe,z) with 

.ft.)-^ ^^^^ . (22, 

Due to the I/2 singularity in the last equation leads to a (pseudo) stationary solution p s . 
which reads after Fourier transformation 

v> f pq(x') ( I x-x' \ \ Id 
Ps{x) = ^dJ D^fM ex H VA" (23) 
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Notice that p s is independent of v. This stationary solution is in accordance with 
correction to this stationary solution reads after Fourier transformation 

= (4^^3/2 / d 3 x'p (x')exp(-(x - x') 2 /4Dt)) 

87T.D y i f-f i \ V £ 



13. The 







x - 




fv 


/ Ai- 



- e-l^l/ferfc (Va7- |X ^ ' j j . (24) 

With other word we find a transient behavior in an intermediate time interval, whereas in 
the long time limit the concentration decays to zero. Such a behavior is similar to the a - 
relaxation in undercooled liquids. 

For d = 2 the hypergeometric function can be expressed by an elementary function. From 
Eq. (jTSjl one gets 

1+ . £ 
p{k,z) = p {k)- 



z + Dk 2 + \ 

Similar to the three dimensional case one has to distinguish between the limiting cases (i) 
kv/z ^> 1 and (ii) kv/z <C 1. Whereas in the first case the initial and the final state 
p p i + pj is reflected, in the second case the solution can be approximately written as 

— * 

p(k, z) = ^ + V(k, z). 

As before for v — > the system tends to a stationary solution of the form 

Psix) = ^ J d 2 x'p Q (x')K (^^) . (25) 

We omit the further details. Instead of that let us consider the (i-dimensional case. Because 
of the factor 1/z in front of the hypergeometric function in Eq. (JT5j) we expect the occurrence 
of a pseudo-stationary for arbitrary d. For t 3> [Dk 2 + A] -1 Eq. (JT5j) can be solved in leading 
order as 

Jd-2 (fcf t) 

p(k, t) ^ Po(k)XA {kvt){d ^ /2[Dk2 + x] . (26) 

This solution can be transformed in the coordinate representation p(x,t). The result is for 
< vt < r =| x — x' I 

= i^pP^ / A ' ( | g- ? ^/. ^W^yfr/Q ■ < 27 > 
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In the opposite case < r < v t we find 

*(*.«) = J(l , t ,^U / A ' ( | e - ffw-^ V fr/a^we ■ < 28 > 

Here /^(x) and K v [x) are Bessel functions. In d = 1 one gets for £ ^> [DA; 2 + A] -1 the 
stationary solution discussed in Eq. (|18|). 

V. CONCLUSIONS 

In this paper we have extended the conventional diffusion reaction system by includ- 
ing non-Markovian memory terms. In particular the growth rate is modified accordingly. 
Namely the particles are permanently subjected to signals of the environment, which acts 
like a background radiation. In case a particles is annihilated a recreation process sets in 
which is triggered by this background signals from the past. Alternatively the model can be 
viewed as a system of traps, distributed randomly at the initial time t = 0. The particles 
immersed in the traps will be activated and released from the traps only when a parti- 
cle is annihilated. Due to the competition between both processes one observes different 
time regimes with a different behavior. In particular the system evolves a transient regime 
with a very slowly varying concentration field which tends under special circumstances to a 
stationary pattern solution in according to previous studies. Such patterns are absent for 
conventional diffusion in an infinite domain without boundary conditions. 
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